Exploration

Load Images

#import imageio.v2 as imageio
import imageio
import pydicom
from pydicom.data import get_testdata_file
import numpy as np

# load a single DICOM image from the scan volume

ct_file = get_testdata_file("CT_small.dcm")
mr_file = get_testdata_file("MR_small.dcm")
im_ct = imageio.imread(ct_file, progress=True) 
im_mr = imageio.imread(mr_file, progress=True)

C:\Users\nasta\AppData\Local\Temp/ipykernel_25900/3685182357.py:11: DeprecationWarning: Starting with ImageIO v3 the behavior of this function will switch to that of iio.v3.imread. To keep the current behavior (and make this warning disappear) use `import imageio.v2 as imageio` or call `imageio.v2.imread` directly.
  im_ct = imageio.imread(ct_file, progress=True)
C:\Users\nasta\AppData\Local\Temp/ipykernel_25900/3685182357.py:12: DeprecationWarning: Starting with ImageIO v3 the behavior of this function will switch to that of iio.v3.imread. To keep the current behavior (and make this warning disappear) use `import imageio.v2 as imageio` or call `imageio.v2.imread` directly.
  im_mr = imageio.imread(mr_file, progress=True)

Medical Image Metadata

Comprising : Patient Demographics and Acquisition Information

ds_ct = pydicom.dcmread(ct_file)
ds_ct

Dataset.file_meta -------------------------------
(0002, 0000) File Meta Information Group Length  UL: 192
(0002, 0001) File Meta Information Version       OB: b'\x00\x01'
(0002, 0002) Media Storage SOP Class UID         UI: CT Image Storage
(0002, 0003) Media Storage SOP Instance UID      UI: 1.3.6.1.4.1.5962.1.1.1.1.1.20040119072730.12322
(0002, 0010) Transfer Syntax UID                 UI: Explicit VR Little Endian
(0002, 0012) Implementation Class UID            UI: 1.3.6.1.4.1.5962.2
(0002, 0013) Implementation Version Name         SH: 'DCTOOL100'
(0002, 0016) Source Application Entity Title     AE: 'CLUNIE1'
-------------------------------------------------
(0008, 0005) Specific Character Set              CS: 'ISO_IR 100'
(0008, 0008) Image Type                          CS: ['ORIGINAL', 'PRIMARY', 'AXIAL']
(0008, 0012) Instance Creation Date              DA: '20040119'
(0008, 0013) Instance Creation Time              TM: '072731'
(0008, 0014) Instance Creator UID                UI: 1.3.6.1.4.1.5962.3
(0008, 0016) SOP Class UID                       UI: CT Image Storage
(0008, 0018) SOP Instance UID                    UI: 1.3.6.1.4.1.5962.1.1.1.1.1.20040119072730.12322
(0008, 0020) Study Date                          DA: '20040119'
(0008, 0021) Series Date                         DA: '19970430'
(0008, 0022) Acquisition Date                    DA: '19970430'
(0008, 0023) Content Date                        DA: '19970430'
(0008, 0030) Study Time                          TM: '072730'
(0008, 0031) Series Time                         TM: '112749'
(0008, 0032) Acquisition Time                    TM: '112936'
(0008, 0033) Content Time                        TM: '113008'
(0008, 0050) Accession Number                    SH: ''
(0008, 0060) Modality                            CS: 'CT'
(0008, 0070) Manufacturer                        LO: 'GE MEDICAL SYSTEMS'
(0008, 0080) Institution Name                    LO: 'JFK IMAGING CENTER'
(0008, 0090) Referring Physician's Name          PN: ''
(0008, 0201) Timezone Offset From UTC            SH: '-0500'
(0008, 1010) Station Name                        SH: 'CT01_OC0'
(0008, 1030) Study Description                   LO: 'e+1'
(0008, 1090) Manufacturer's Model Name           LO: 'RHAPSODE'
(0009, 0010) Private Creator                     LO: 'GEMS_IDEN_01'
(0009, 1001) [Full fidelity]                     LO: 'GE_GENESIS_FF'
(0009, 1002) [Suite id]                          SH: 'CT01'
(0009, 1004) [Product id]                        SH: 'HiSpeed CT/i'
(0009, 1027) [Image actual date]                 SL: 862399669
(0009, 1030) [Service id]                        SH: ''
(0009, 1031) [Mobile location number]            SH: ''
(0009, 10e6) [Genesis Version - now]             SH: '05'
(0009, 10e7) [Exam Record checksum]              UL: 973283917
(0009, 10e9) [Actual series data time stamp]     SL: 862399669
(0010, 0010) Patient's Name                      PN: 'CompressedSamples^CT1'
(0010, 0020) Patient ID                          LO: '1CT1'
(0010, 0030) Patient's Birth Date                DA: ''
(0010, 0040) Patient's Sex                       CS: 'O'
(0010, 1002)  Other Patient IDs Sequence  2 item(s) ---- 
   (0010, 0020) Patient ID                          LO: 'ABCD1234'
   (0010, 0022) Type of Patient ID                  CS: 'TEXT'
   ---------
   (0010, 0020) Patient ID                          LO: '1234ABCD'
   (0010, 0022) Type of Patient ID                  CS: 'TEXT'
   ---------
(0010, 1010) Patient's Age                       AS: '000Y'
(0010, 1030) Patient's Weight                    DS: '0.0'
(0010, 21b0) Additional Patient History          LT: ''
(0011, 0010) Private Creator                     LO: 'GEMS_PATI_01'
(0011, 1010) [Patient Status]                    SS: 0
(0018, 0010) Contrast/Bolus Agent                LO: 'ISOVUE300/100'
(0018, 0022) Scan Options                        CS: 'HELICAL MODE'
(0018, 0050) Slice Thickness                     DS: '5.0'
(0018, 0060) KVP                                 DS: '120.0'
(0018, 0088) Spacing Between Slices              DS: '5.0'
(0018, 0090) Data Collection Diameter            DS: '480.0'
(0018, 1020) Software Versions                   LO: '05'
(0018, 1040) Contrast/Bolus Route                LO: 'IV'
(0018, 1100) Reconstruction Diameter             DS: '338.6716'
(0018, 1110) Distance Source to Detector         DS: '1099.3100585938'
(0018, 1111) Distance Source to Patient          DS: '630.0'
(0018, 1120) Gantry/Detector Tilt                DS: '0.0'
(0018, 1130) Table Height                        DS: '133.699997'
(0018, 1150) Exposure Time                       IS: '1601'
(0018, 1151) X-Ray Tube Current                  IS: '170'
(0018, 1152) Exposure                            IS: '170'
(0018, 1160) Filter Type                         SH: 'LARGE BOWTIE FIL'
(0018, 1190) Focal Spot(s)                       DS: '0.7'
(0018, 1210) Convolution Kernel                  SH: 'STANDARD'
(0018, 5100) Patient Position                    CS: 'FFS'
(0019, 0010) Private Creator                     LO: 'GEMS_ACQU_01'
(0019, 1002) [Detector Channel]                  SL: 912
(0019, 1003) [Cell number at Theta]              DS: '373.75'
(0019, 1004) [Cell spacing]                      DS: '1.0166'
(0019, 100f) [Horiz. Frame of ref.]              DS: '955.799988'
(0019, 1011) [Series contrast]                   SS: 2
(0019, 1013) [Start number for baseline]         SS: 0
(0019, 1014) [End number for baseline]           SS: 0
(0019, 1015) [Start number for enhanced scans]   SS: 0
(0019, 1016) [End number for enhanced scans]     SS: 0
(0019, 1017) [Series plane]                      SS: 2
(0019, 1018) [First scan ras]                    LO: 'S'
(0019, 1019) [First scan location]               DS: '7.79187'
(0019, 101a) [Last scan ras]                     LO: 'I'
(0019, 101b) [Last scan loc]                     DS: '-320.197968'
(0019, 101e) [Display field of view]             DS: '0.0'
(0019, 1023) [Table Speed [mm/rotation]]         DS: '5.0'
(0019, 1024) [Mid Scan Time [sec]]               DS: '17.784578'
(0019, 1025) [Mid scan flag]                     SS: 1
(0019, 1026) [Tube Azimuth [degree]]             SL: 0
(0019, 1027) [Rotation Speed [msec]]             DS: '1.0'
(0019, 102a) [x-ray On position]                 DS: '178.079926'
(0019, 102b) [x-ray Off position]                DS: '3994.299316'
(0019, 102c) [Number of triggers]                SL: 10431
(0019, 102e) [Angle of first view]               DS: '-718.079956'
(0019, 102f) [Trigger frequency]                 DS: '984.0'
(0019, 1039) [SFOV Type]                         SS: 16
(0019, 1040) [Stat recon flag]                   SS: 0
(0019, 1041) [Compute type]                      SS: 1
(0019, 1042) [Segment Number]                    SS: 0
(0019, 1043) [Total Segments Required]           SS: 0
(0019, 1044) [Interscan delay]                   DS: '1.0'
(0019, 1047) [View compression factor]           SS: 1
(0019, 104a) [Total no. of ref channels]         SS: 6
(0019, 104b) [Data size for scan data]           SL: 10431
(0019, 1052) [Recon post proc. Flag]             SS: 1
(0019, 1057) [CT water number]                   SS: -95
(0019, 1058) [CT bone number]                    SS: 0
(0019, 105e) [Number of channels (1...512)]      SL: 763
(0019, 105f) [Increment between channels]        SL: 1
(0019, 1060) [Starting view]                     SL: 1969
(0019, 1061) [Number of views]                   SL: 1576
(0019, 1062) [Increment between views]           SL: 1
(0019, 106a) [Dependent on #views processed]     SS: 4
(0019, 106b) [Field of view in detector cells]   SS: 852
(0019, 1070) [Value of back projection button]   SS: 1
(0019, 1071) [Set if fatq estimates were used]   SS: 0
(0019, 1072) [Z chan avg over views]             DS: '0.0'
(0019, 1073) [Avg of left ref chans over views]  DS: '0.0'
(0019, 1074) [Max left chan over views]          DS: '0.0'
(0019, 1075) [Avg of right ref chans over views] DS: '0.0'
(0019, 1076) [Max right chan over views]         DS: '0.0'
(0019, 10da) [Reference channel used]            SS: 0
(0019, 10db) [Back projector coefficient]        DS: '0.0'
(0019, 10dc) [Primary speed correction used]     SS: 1
(0019, 10dd) [Overrange correction used]         SS: 1
(0019, 10de) [Dynamic Z alpha value]             DS: '0.0'
(0020, 000d) Study Instance UID                  UI: 1.3.6.1.4.1.5962.1.2.1.20040119072730.12322
(0020, 000e) Series Instance UID                 UI: 1.3.6.1.4.1.5962.1.3.1.1.20040119072730.12322
(0020, 0010) Study ID                            SH: '1CT1'
(0020, 0011) Series Number                       IS: '1'
(0020, 0012) Acquisition Number                  IS: '2'
(0020, 0013) Instance Number                     IS: '1'
(0020, 0032) Image Position (Patient)            DS: [-158.135803, -179.035797, -75.699997]
(0020, 0037) Image Orientation (Patient)         DS: [1.000000, 0.000000, 0.000000, 0.000000, 1.000000, 0.000000]
(0020, 0052) Frame of Reference UID              UI: 1.3.6.1.4.1.5962.1.4.1.1.20040119072730.12322
(0020, 0060) Laterality                          CS: ''
(0020, 1040) Position Reference Indicator        LO: 'SN'
(0020, 1041) Slice Location                      DS: '-77.2040634155'
(0020, 4000) Image Comments                      LT: 'Uncompressed'
(0021, 0010) Private Creator                     LO: 'GEMS_RELA_01'
(0021, 1003) [Series from which Prescribed]      SS: 2
(0021, 1005) [Genesis Version - now]             SH: '05'
(0021, 1007) [Series Record checksum]            UL: 1605775145
(0021, 1015) [Unknown]                           US: 24078
(0021, 1016) [Unknown]                           SS: 2
(0021, 1018) [Genesis version - Now]             SH: '05'
(0021, 1019) [Acq recon record checksum]         UL: 750675506
(0021, 1037) [Screen Format]                     SS: 16
(0021, 104a) [Anatomical reference for scout]    LO: ''
(0021, 1090) [Tube focal spot position]          SS: 7400
(0021, 1091) [Biopsy position]                   SS: 0
(0021, 1092) [Biopsy T location]                 FL: 0.0
(0021, 1093) [Biopsy ref location]               FL: 0.0
(0023, 0010) Private Creator                     LO: 'GEMS_STDY_01'
(0023, 1070) [Start time(secs) in first axial]   FD: 862399761.111079
(0023, 1074) [No. of updates to header]          SL: 1
(0023, 107d) [Indicates study has complete info  SS: 0
(0025, 0010) Private Creator                     LO: 'GEMS_SERS_01'
(0025, 1006) [Last pulse sequence used]          SS: 0
(0025, 1007) [Images in Series]                  SL: 44
(0025, 1010) [Landmark Counter]                  SL: 0
(0025, 1011) [Number of Acquisitions]            SS: 0
(0025, 1017) [Series Complete Flag]              SL: 0
(0025, 1018) [Number of images archived]         SL: 0
(0025, 1019) [Last image number used]            SL: 4
(0025, 101a) [Primary Receiver Suite and Host]   SH: ''
(0027, 0010) Private Creator                     LO: 'GEMS_IMAG_01'
(0027, 1006) [Image archive flag]                SL: 1
(0027, 1010) [Scout Type]                        SS: 0
(0027, 101c) [Vma mamp]                          SL: 150
(0027, 101d) [Vma phase]                         SS: 1
(0027, 101e) [Vma mod]                           SL: 24
(0027, 101f) [Vma clip]                          SL: 129
(0027, 1020) [Smart scan ON/OFF flag]            SS: 1
(0027, 1030) [Foreign Image Revision]            SH: ''
(0027, 1035) [Plane Type]                        SS: 2
(0027, 1040) [RAS letter of image location]      SH: 'I'
(0027, 1041) [Image location]                    FL: -77.20406341552734
(0027, 1042) [Center R coord of plane image]     FL: -11.199999809265137
(0027, 1043) [Center A coord of plane image]     FL: 9.699999809265137
(0027, 1044) [Center S coord of plane image]     FL: -75.69999694824219
(0027, 1045) [Normal R coord]                    FL: 0.0
(0027, 1046) [Normal A coord]                    FL: 0.0
(0027, 1047) [Normal S coord]                    FL: -1.0
(0027, 1048) [R Coord of Top Right Corner]       FL: -180.53579711914062
(0027, 1049) [A Coord of Top Right Corner]       FL: 179.03579711914062
(0027, 104a) [S Coord of Top Right Corner]       FL: -75.69999694824219
(0027, 104b) [R Coord of Bottom Right Corner]    FL: -180.53579711914062
(0027, 104c) [A Coord of Bottom Right Corner]    FL: -159.63580322265625
(0027, 104d) [S Coord of Bottom Right Corner]    FL: -75.69999694824219
(0027, 1050) [Scan Start Location]               FL: -63.19999694824219
(0027, 1051) [Scan End Location]                 FL: -116.20304870605469
(0027, 1052) [RAS letter for side of image]      SH: 'L'
(0027, 1053) [RAS letter for anterior/posterior] SH: 'A'
(0027, 1054) [RAS letter for scout start loc]    SH: 'I'
(0027, 1055) [RAS letter for scout end loc]      SH: 'I'
(0028, 0002) Samples per Pixel                   US: 1
(0028, 0004) Photometric Interpretation          CS: 'MONOCHROME2'
(0028, 0010) Rows                                US: 128
(0028, 0011) Columns                             US: 128
(0028, 0030) Pixel Spacing                       DS: [0.661468, 0.661468]
(0028, 0100) Bits Allocated                      US: 16
(0028, 0101) Bits Stored                         US: 16
(0028, 0102) High Bit                            US: 15
(0028, 0103) Pixel Representation                US: 1
(0028, 0120) Pixel Padding Value                 SS: -2000
(0028, 1052) Rescale Intercept                   DS: '-1024.0'
(0028, 1053) Rescale Slope                       DS: '1.0'
(0029, 0010) Private Creator                     LO: 'GEMS_IMPS_01'
(0029, 1004) [Lower range of Pixels1]            SL: 0
(0029, 1005) [Lower range of Pixels1]            DS: '0.0'
(0029, 1006) [Lower range of Pixels1]            DS: '0.0'
(0029, 1007) [Lower range of Pixels1]            SL: 87
(0029, 1008) [Lower range of Pixels1]            SH: ''
(0029, 1009) [Lower range of Pixels1]            SH: ''
(0029, 100a) [Lower range of Pixels1]            SS: 764
(0029, 1026) [Version of the hdr struct]         SS: 2
(0029, 1034) [Advantage comp. Overflow]          SL: 0
(0029, 1035) [Advantage comp. Underflow]         SL: 0
(0043, 0010) Private Creator                     LO: 'GEMS_PARM_01'
(0043, 1010) [Window value]                      US: 400
(0043, 1011) [Total input views]                 US: 10431
(0043, 1012) [X-ray chain]                       SS: [14, 2, 3]
(0043, 1013) [Decon kernel parameters]           SS: [107, 21, 4, 2, 20]
(0043, 1014) [Calibration parameters]            SS: [4, 4, 5]
(0043, 1015) [Total output views]                SS: 10431
(0043, 1016) [Number of overranges]              SS: 0
(0043, 1017) [IBH image scale factors]           DS: '0.095'
(0043, 1018) [BBH coefficients]                  DS: [0.085000, 1.102000, 0.095000]
(0043, 1019) [Number of BBH chains to blend]     SS: 350
(0043, 101a) [Starting channel number]           SL: 7
(0043, 101b) [Ppscan parameters]                 SS: [0, 0, 0, 0, 0]
(0043, 101c) [GE image integrity]                SS: 0
(0043, 101d) [Level value]                       SS: 40
(0043, 101e) [Delta Start Time [msec]]           DS: '2.0'
(0043, 101f) [Max overranges in a view]          SL: 0
(0043, 1020) [Avg overranges all views]          DS: '0.0'
(0043, 1021) [Corrected after glow terms]        SS: 0
(0043, 1025) [Reference channels]                SS: [1, 2, 3, 748, 749, 750]
(0043, 1026) [No views ref chans blocked]        US: [0, 1, 1, 0, 0, 0]
(0043, 1027) [Scan Pitch Ratio]                  SH: '/1.0:1'
(0043, 1028) [Unique image iden]                 OB: Array of 80 elements
(0043, 1029) [Histogram tables]                  OB: Array of 2068 elements
(0043, 102a) [User defined data]                 OB: Array of 40 elements
(0043, 102b) [Private Scan Options]              SS: [4, 4, 0, 0]
(0043, 1031) [Recon Center Coordinates]          DS: [-11.200000, 9.700000]
(0043, 1040) [Trigger on position]               FL: 178.07992553710938
(0043, 1041) [Degree of rotation]                FL: 3816.219482421875
(0043, 1042) [DAS trigger source]                SL: 0
(0043, 1043) [DAS fpa gain]                      SL: 0
(0043, 1044) [DAS output source]                 SL: 1
(0043, 1045) [DAS ad input]                      SL: 1
(0043, 1046) [DAS cal mode]                      SL: 3
(0043, 1047) [DAS cal frequency]                 SL: -1
(0043, 1048) [DAS reg xm]                        SL: 1
(0043, 1049) [DAS auto zero]                     SL: 0
(0043, 104a) [Starting channel of view]          SS: 1
(0043, 104b) [DAS xm pattern]                    SL: 0
(0043, 104c) [TGGC trigger mode]                 SS: 0
(0043, 104d) [Start scan to X-ray on delay]      FL: 0.0
(0043, 104e) [Duration of X-ray on]              FL: 10.60060977935791
(7fe0, 0010) Pixel Data                          OW: Array of 32768 elements
(fffc, fffc) Data Set Trailing Padding           OB: Array of 126 elements

Plot Images

Matplotlib's imshow() function </b> cmap='gray' means gray color mappings for each value

import matplotlib.pyplot as plt

# Draw the image in grayscale
plt.imshow(im_ct, cmap='gray')

# Render the image
plt.show()

Load Volumes

ImageIO's volread() function
Load multi-dimensional datasets and create 3D volumes from a folder of images

# Path to the folder containing the DICOM files
folder_path = 'Data/Brain/SE000001/MR000000'

# Load the brain data folder
vol = imageio.volread(folder_path)

Reading DICOM (examining files): 1/27 files (3.7%27/27 files (100.0%)
  Found 1 correct series.
Reading DICOM (loading data): 27/27  (100.0%)

Slice 3D Images

Slicing 3D and 4D images into many 2D frames

# Plot the first slice of the volume (2D frame)
plt.imshow(vol[0, :, :], cmap='gray') 

# Render the image
plt.show()

Plot Other Views

Slicing along different axes

# Compute aspect ratios
d0, d1, d2 = vol.meta['sampling']
asp1 = d0 / d2
# Plot the other view slice of the volume
plt.imshow(vol[:, 128, :], cmap='gray',aspect=asp1)
# Render the image
plt.show()

Image Comparison

Translation

Shifting or moving an image from one place to another within a picture.

import scipy.ndimage as ndi

#formatting method for the plots in this notebook:
def format_and_render_plot():
    '''
    Custom function to simplify common formatting operations for exercises. Operations include: 
    1. Turning off axis grids.
    2. Calling `plt.tight_layout` to improve subplot spacing.
    3. Calling `plt.show()` to render plot.
    '''
    fig = plt.gcf()
    for ax in fig.axes:
        ax.axis('off')    
    plt.tight_layout()
    plt.show()

# save the middle slice as separat image
middle_slice = vol.shape[0] // 2            # // is floor division
im = vol[middle_slice,:,:]

# Translate the brain towards the center
xfm = ndi.shift(im, shift=(10, 15))

# Plot the original and adjusted images
fig, axes = plt.subplots(nrows=1, ncols=2)
axes[0].imshow(im, cmap='gray')
axes[0].set_title('original', fontweight ="bold")
axes[1].imshow(xfm, cmap='gray')
axes[1].set_title('translated', fontweight ="bold")
format_and_render_plot()

Rotation

Rotating an image from its center by the specified degrees from the right horizontal axis.

# Rotate the shifted image
xfm = ndi.rotate(xfm, angle=-30, reshape=False)

# Plot the original and transformed images
fig, axes = plt.subplots(nrows=1, ncols=2)
axes[0].imshow(im, cmap='gray')
axes[0].set_title('original', fontweight ="bold")
axes[1].imshow(xfm, cmap='gray')
axes[1].set_title('rotated', fontweight ="bold")
format_and_render_plot()

Affine Transformation

Scaling (changing size), rotation, translation (shifting), and shearing (changing angles) using matrix multiplication

# Define the affine transform matrix
mat = np.array([[0.8, -0.4, 90], [0.4, 0.8, -6.0], [0, 0, 1]])
print(mat)

# Apply the affine transform matrix to image data
xfm = ndi.affine_transform(im, mat)

# Plot the original and transformed images
fig, axes = plt.subplots(nrows=1, ncols=2)
axes[0].imshow(im, cmap='gray')
axes[0].set_title('original', fontweight ="bold")
axes[1].imshow(xfm, cmap='gray')
axes[1].set_title('transformed', fontweight ="bold")
format_and_render_plot()

[[ 0.8 -0.4 90. ]
 [ 0.4  0.8 -6. ]
 [ 0.   0.   1. ]]

Resampling

• Downsampling: combining pixel data to decrease size
• Upsampling: distributing pixel data to increase size

# Resample image
im_dn = ndi.zoom(im, zoom=0.25)
im_up = ndi.zoom(im, zoom=4.00)

# Plot the images
fig, axes = plt.subplots(1, 2)
axes[0].imshow(im_dn, cmap='gray')
axes[0].set_title('downsampled', fontweight ="bold")
axes[1].imshow(im_up, cmap='gray')
axes[1].set_title('upsampled', fontweight ="bold")
format_and_render_plot()

Interpolation

Estimating new pixel intensities after an image transformation

# Upsample "im" by a factor of 4
up0 = ndi.zoom(im, zoom=512/128, order=0)
up5 = ndi.zoom(im, zoom=512/128, order=5)

# Print original and new shape
print('Original shape:', im.shape)
print('Upsampled shape:', up5.shape)

# Plot close-ups of the new images
fig, axes = plt.subplots(1, 2)
axes[0].imshow(up0[128:256, 128:256], cmap='gray')
axes[0].set_title('no interpolation', fontweight ="bold")
axes[1].imshow(up5[128:256, 128:256], cmap='gray')
axes[1].set_title('with interpolation', fontweight ="bold")
format_and_render_plot()

Original shape: (256, 256)
Upsampled shape: (1024, 1024)

Mean Absolute Error (MAE)

Mean intensity differences between two images, with higher values indicating greater divergence

# Load the dcm file (image)
im1 = im

# Apply the affine transform matrix to image data
xfm = ndi.shift(im, shift=(-12, -8))
im2 = xfm

# Calculate image difference
err = im1 - im2

# Plot the difference
fig, axes = plt.subplots(1, 3)
axes[0].imshow(im1, cmap='gray')
axes[0].set_title('original', fontweight ="bold")
axes[1].imshow(im2, cmap='gray')
axes[1].set_title('shifted', fontweight ="bold")
axes[2].imshow(err, cmap='seismic', vmin=-200, vmax=200)
axes[2].set_title('error map', fontweight ="bold")
format_and_render_plot()

# Calculate absolute image difference
abs_err = np.absolute(im1 - im2)

# Plot the difference
fig, axes = plt.subplots(1, 3)
axes[0].imshow(im1, cmap='gray')
axes[0].set_title('original', fontweight ="bold")
axes[1].imshow(im2, cmap='gray')
axes[1].set_title('shifted', fontweight ="bold")
axes[2].imshow(abs_err, cmap='seismic', vmin=-200, vmax=200)
axes[2].set_title('error map (absolute)', fontweight ="bold")
format_and_render_plot()

# Calculate mean absolute error
mean_abs_err = np.mean(np.abs(im1 - im2))
print('MAE:', mean_abs_err)

MAE: 25730.311309814453

Structural Similarity (SSIM)

Structural similarity between two images perceived changes in:
structural information, luminance, and contrast that humans typically notice.

from skimage.metrics import structural_similarity as ssim

# Compute SSIM between im and xfm
data_range=xfm.max() - xfm.min()
ssim_index = ssim(im, xfm, data_range=data_range)

print(f'SSIM: {ssim_index}')

SSIM: 0.47951331548119075

Intersection of The Union (IoU)

The number of pixels filled in both images (the intersection) out of the number of pixels filled in either image (the union).

def intersection_of_union(im1, im2):
    i = np.logical_and(im1, im2)
    u = np.logical_or(im1, im2)
    return i.sum() / u.sum()

# Try some other paramters by yourself
xfm = ndi.shift(im, shift=(-10, -10))
xfm = ndi.rotate(xfm, angle=-15, reshape=False)
intersection_of_union(xfm, im)

0.8135092272202998

Combining Two Images

The number of pixels filled in both images (the intersection) out of the number of pixels filled in either image (the union).

import SimpleITK as sitk
import Utilities.gui
img1 = sitk.ReadImage('../ImageAnalysis/Data/Brain/training_001_mr_T1.mha') #MRI image
img2_original = sitk.ReadImage('../ImageAnalysis/Data/Brain/training_001_ct.mha')  #CT image
img2 = sitk.Resample(img2_original, img1)

# Obtain foreground masks for the two images using Otsu thresholding, we use these later on.
msk1 = sitk.OtsuThreshold(img1,0,1)
msk2 = sitk.OtsuThreshold(img2,0,1)

Utilities.gui.MultiImageDisplay(image_list = [img1, img2],                   
                      title_list = ['image1', 'image2'],
                      figure_size=(9,3))

<Utilities.gui.MultiImageDisplay at 0x2a5c2f58280>

Masks and Filters

Load Images

#  now we load test data files (CT and MR)
ct_file = get_testdata_file("CT_small.dcm")

# Load the CT file (image)
im = imageio.imread(ct_file)

C:\Users\nasta\AppData\Local\Temp/ipykernel_25900/2938429760.py:5: DeprecationWarning: Starting with ImageIO v3 the behavior of this function will switch to that of iio.v3.imread. To keep the current behavior (and make this warning disappear) use `import imageio.v2 as imageio` or call `imageio.v2.imread` directly.
  im = imageio.imread(ct_file)

Intensity

Intensity Normalization

# min-max normalisation
im_old = im                                     # save original image

print('Min. value before normalization:', im_old.min())
print('Max value before normalization:', im_old.max())

im = (im - im.min()) / (im.max() - im.min())    # normalise to 0-1 range

print('Min. value after normalization:', im.min())
print('Max value after normalization:', im.max())

Min. value before normalization: -896
Max value before normalization: 1167
Min. value after normalization: 0.0
Max value after normalization: 1.0

def format_and_render_plot(axis=False, legend=False):
    '''
    Custom function to simplify common formatting operations for exercises. Operations include: 
    1. Turning off axis grids and legends, if not explicitly requested.
    2. Calling `plt.tight_layout` to improve subplot spacing.
    3. Calling `plt.show()` to render plot.
    '''
    fig = plt.gcf()
    for ax in fig.axes:
        if not axis:
            ax.axis('off')
        if legend:  
            ax.legend(loc='center right')  
    plt.tight_layout()
    plt.show()
fig, axes = plt.subplots(1, 2, sharex=True)
axes[0].imshow(im_old, cmap='gray')
axes[0].set_title('without normalization', fontweight ="bold")
axes[1].imshow(im, cmap='gray')
axes[1].set_title('with normalization', fontweight ="bold")
format_and_render_plot()

Histogram

The distribution of values in your image by binning each element by its intensity then measuring the size of each bin.

CDF (Cumulative distribution function)

The frequency with which a given range of pixel intensities occurs.

# Create a 256-bin histogram, binned at each possible value
hist = ndi.histogram(im, min=im.min(), max=im.max(), bins=256)

# Create a cumulative distribution function
cdf = hist.cumsum() / hist.sum()

# Plot the histogram and CDF
fig, axes = plt.subplots(2, 1, sharex=True)         # sharex=True shares the x-axis between the top and bottom subplot
axes[0].plot(hist, label='Histogram')               # label='Histogram' labels this line as "Histogram" for the legend
axes[1].plot(cdf, label='CDF')                      # label='CDF' labels this line as "CDF" for the legend
format_and_render_plot(axis=True, legend=True)      # axis=True turns on axis grids for the plot; legend=True turns on the legend

Create a Mask

Binary arrays for removing or selecting specific parts of an image by applying one or more logical operations to an image.

# Create soft tissue (st) and bone masks
# Try different intervals by yourself
mask_st = (im >= 0.2) & (im < 0.52)     # soft tissue mask
mask_bone = im >= 0.53                  # bone mask

# Plot the skin (0) and bone (1) masks
fig, axes = plt.subplots(1,2)
axes[0].imshow(mask_st, cmap='gray')
axes[0].set_title('soft tissue map', fontweight ="bold")
axes[1].imshow(mask_bone, cmap='gray')
axes[1].set_title('bone map', fontweight ="bold")
format_and_render_plot()

Apply a Mask

Although masks are binary, they can be applied to images to filter out pixels where the mask is False.

NumPy's where() function is a flexible way of applying masks. It takes three arguments:

np.where(condition, x, y)

condition, x and y can be either arrays or single values. This allows you to pass through original image values while setting masked values to 0.

# Screen out non-bone pixels from "im"
im_bone = np.where(mask_bone, im, 0)  #extract pixels from im which locate in mask_bone!

# Plot masked image and histogram
plt.imshow(im_bone, cmap='gray', label='Bone pixels')

<matplotlib.image.AxesImage at 0x2a5c6fd0520>

Tune a Mask

Imperfect masks can be tuned through the addition and subtraction of pixels. example methods for accomplishing these ends:

• binary_dilation: Add pixels along edges
• binary_closing: Dilate then erode, "filling in" holes

# Create and tune bone mask
mask_dilate = ndi.binary_dilation(mask_bone, iterations=5)
mask_closed = ndi.binary_closing(mask_bone, iterations=5)

# Plot masked images
fig, axes = plt.subplots(1,3)
axes[0].imshow(mask_bone)
axes[0].set_title('original', fontweight ="bold")
axes[1].imshow(mask_dilate)
axes[1].set_title('dilated', fontweight ="bold")
axes[2].imshow(mask_closed)
axes[2].set_title('closed', fontweight ="bold")
format_and_render_plot()

Filter Convolutions

Transforming images based on intensity values surrounding a pixel

# Set filter weights
weights = [[0.11, 0.11, 0.11],
           [0.11, 0.11, 0.11], 
           [0.11, 0.11, 0.11]]

# Convolve the image with the filter
im_filt = ndi.convolve(im, weights)

# Plot the images
fig, axes = plt.subplots(1,2)
axes[0].imshow(im, cmap='gray')
axes[0].set_title('original', fontweight ="bold")
axes[1].imshow(im_filt, cmap='gray')
axes[1].set_title('with filter', fontweight ="bold")
format_and_render_plot()

Smoothing

Improve the signal-to-noise ratio of an image by blurring out small variations in intensity. The Gaussian filter is excellent for this: it is a circular (or spherical) smoothing kernel that weights nearby pixels higher than distant ones.

# Smooth "im" with Gaussian filters
# Try different sigmas  by yourself
im_s1 = ndi.gaussian_filter(im, sigma=1)
im_s3 = ndi.gaussian_filter(im, sigma=3)

# Draw bone masks of each image
fig, axes = plt.subplots(1,3)
axes[0].imshow(im, cmap='gray')
axes[0].set_title('original', fontweight ="bold")
axes[1].imshow(im_s1, cmap='gray')
axes[1].set_title('Gaussian (sigma=1)', fontweight ="bold")
axes[2].imshow(im_s3, cmap='gray')
axes[2].set_title('Gaussian (sigma=3)', fontweight ="bold")
format_and_render_plot()

Detect Edges

Applying filters which pattern is a change in intensity along a plane.

# Set weights to detect vertical edges
weights = [[1, 0, -1],
           [1, 0, -1],
           [1, 0, -1]]

# Convolve "im" with filter weights
edges = ndi.convolve(im, weights)

# Draw the image in color
fig, axes = plt.subplots(1,2)
axes[0].imshow(im, cmap='gray')
axes[0].set_title('original', fontweight ="bold")
axes[1].imshow(edges, cmap='seismic', vmin=-0.8, vmax=0.8)
axes[1].set_title('edge filter', fontweight ="bold")
format_and_render_plot()

Measurement

Load Image/Volume

# Load the directory (volume)
folder_path = 'Data/Brain/SE000001/MR000000'

# Load the volume
vol = imageio.volread(folder_path)
print(vol.shape)

# save the middle slice as separat image
middle_slice = vol.shape[0] // 2            # // is floor division
im = vol[middle_slice,:,:]
# min-max normalisation
im_old = im                                     # save original image for later
im = (im - im.min()) / (im.max() - im.min())    # normalise the image, range 0 - 1

Reading DICOM (examining files): 1/27 files (3.7%27/27 files (100.0%)
  Found 1 correct series.
Reading DICOM (loading data): 27/27  (100.0%)
(27, 256, 256)

Segmentation

Apply a mask and fill any gaps in the mask

# Smooth intensity values
im_filt = ndi.median_filter(im, size=3)         # size = 3 means 3x3x3 neighbourhood

# Select high-intensity pixels
mask_start = np.where(im_filt > 0.3, 1, 0)      # mask_start is a boolean array
mask = ndi.binary_closing(mask_start)           # fill holes

# Label the objects in "mask"
labels, nlabels = ndi.label(mask)               # labels: each object has a unique number, nlabels: number of objects

# Create a `labels` overlay
overlay = np.where(labels > 0, labels, np.nan)

# Use imshow to plot the overlay
fig, axes = plt.subplots(1, 2)
axes[0].imshow(im, cmap='gray')
axes[0].set_title('original', fontweight ="bold")
axes[1].imshow(im, cmap='gray')                             # show image first
axes[1].imshow(overlay, cmap='rainbow', alpha=0.6)          # show overlay second; alpha controls transparency
axes[1].set_title('with segmentation', fontweight ="bold")
format_and_render_plot()

Select Objects

Pick up whole sets of pixels at a time

# Label the image "mask"
labels, nlabels = ndi.label(mask)

# Select brain label
brain_val = 2                                           # 2 is the label for the brain (see plot above)
brain_mask = np.where(labels == brain_val, 1, np.nan)   # create brain mask

# Overlay selected label
fig, axes = plt.subplots(1, 2)
axes[0].imshow(im, cmap='gray')
axes[0].set_title('original', fontweight ="bold")
axes[1].imshow(im, cmap='gray')
axes[1].imshow(brain_mask, cmap='rainbow', alpha=0.6)
axes[1].set_title('with segmentation', fontweight ="bold")
format_and_render_plot()

Extract Objects

Crop images so that they only include the object of interest.

brain_mask=brain_mask.astype(np.int64)

# Find bounding box
bboxes =ndi.find_objects(brain_mask)                # returns a list of tuples, each tuple has 3 slices
print('Number of objects:', len(bboxes))            # number of objects found
print('Indices for first box:', bboxes[0])          # print indices for first box

# Crop to index 0
im_brain = im[bboxes[0]]                            # crop the original image to the bounding box

# Plot the cropped image
fig, axes = plt.subplots(1, 2)
axes[0].imshow(im, cmap='gray')
axes[0].set_title('original', fontweight ="bold")
axes[1].imshow(im_brain, cmap='gray')
axes[1].set_title('cropped', fontweight ="bold")
format_and_render_plot()

Number of objects: 1
Indices for first box: (slice(50, 234, None), slice(61, 193, None))

Measure Variance

Tailor measurements to specific sets of pixels:

# Variance for all pixels
var_all = ndi.variance(vol, labels=None, index=None)            # labels=None means all pixels, index=None means all objects
print('All pixels:', var_all)

# Variance for labeled pixels
var_labels = ndi.variance(vol, labels, index=None)              # labels=labels means only labeled pixels, index=None means all objects
print('Labeled pixels:', var_labels)

# Variance for each object
var_objects = ndi.variance(vol, labels, index=[1,2])            # labels=labels means only labeled pixels, index=[1,2] means only objects 1 and 2
print('Brain matter:', var_objects[0])
print('Other tissue:', var_objects[1])

All pixels: 15281.189220274951
Labeled pixels: 18188.560688187546
Brain matter: 13466.402775068886
Other tissue: 16319.82340391041

Separate Histograms

Wide distribution of intensity values and more variance of multiple tissue types

# Create histograms for selected pixels
hist1 = ndi.histogram(vol, min=0, max=255, bins=256)
hist2 = ndi.histogram(vol, 0, 255, 256, labels=labels)
hist3 = ndi.histogram(vol, 0, 255, 256, labels=labels, index=1)


# Plot the histogram and CDF
fig, axes = plt.subplots(3, 1, sharex=True)                 # sharex=True shares the x-axis between the top and bottom subplot
axes[0].plot(hist1 / hist1.sum(), label='All pixels')
axes[1].plot(hist2 / hist2.sum(), label='All labeled pixels')
axes[2].plot(hist3 / hist3.sum(), label='Brain matter')
format_and_render_plot(axis=True, legend=True)              # axis=True turns on axis grids for the plot; legend=True turns on the legend

Calculate Distance

the distance from each pixel to a given point

# Calculate distances
brain = np.where(labels == 2, 1, 0)                                             # create brain mask
dists = ndi.distance_transform_edt(brain, sampling=vol.meta['sampling'][1:3])   # calculate distances

# Report on distances
print('Max distance (mm):', ndi.maximum(dists))
print('Max location:', ndi.maximum_position(dists))

# Plot overlay of distances
overlay = np.where(dists > 0, dists, np.nan)                                    # create overlay

fig, axes = plt.subplots(1, 2)
axes[0].imshow(im, cmap='gray')
axes[0].set_title('original', fontweight ="bold")
axes[1].imshow(im, cmap='gray')
axes[1].imshow(overlay, cmap='hot', alpha=0.75)
axes[1].set_title('with distances', fontweight ="bold")
format_and_render_plot()

Max distance (mm): 12.1875
Max location: (133, 175)

Center-Of-Mass (COM)

The coordinates for the point of an object with higher values pulling the center closer to it.

# Extract centers of mass for objects 1 and 2
coms = ndi.center_of_mass(vol, labels, index=[1,2])
print('Label 1 center:', coms[0])
print('Label 2 center:', coms[1])

# Add marks to plot
for c0, c1, c2 in coms:
    plt.scatter(c2, c1, s=100, marker='o')
plt.show()

Label 1 center: (8.856682792208547, 49.44438084201855, 125.03554483621987)
Label 2 center: (11.713445688803832, 119.3292456288401, 129.36906622106693)